library("dplyr")
library("naniar")
library("sjPlot")
library("sjlabelled")
library("sjmisc")
library("ggplot2")
library("broom")
library("mediation")
library("knitr")
library("interplot")
library("readxl")
library("interactions")
library("psych")
library("Hmisc")
library("sensemakr")
library("MASS")
library("stargazer")

##set directory 

setwd("~/Downloads/Business/EduHK/1.COVID project/Jan 12.2022/Codes")

df <- read.csv("covidimg")

##data cleaning

na_strings <- c(96,97,98,99)

df[,c(1,3,4,6)] <- list(NULL)

df2 <- df %>%
  filter(TIMES_3 > 300, TIMES_3 <7200)%>%
  replace_with_na_all(condition = ~.x %in% na_strings)%>%
  mutate_if(is.character,as.numeric)

df2$Anxiety_M <- rowMeans(df2[, c("Q23_1", "Q23_2", "Q23_3", "Q23_4", "Q23_5")], na.rm= TRUE)

df3 <- df2%>%
  mutate(xHQ41 =recode(Q41,`1` = 1, `2` = 1, `3` = 1, `4` = 2, `7` = 3, `6` = 3, `5` = 3,))



##Alpha test

alpha(df3[,c("Q23_1", "Q23_2", "Q23_3", "Q23_4", "Q23_5")])

#Summary 

summary(df3$Anxiety_M)
sd(df3$Anxiety_M, na.rm=T)

#Create interaction terms 

df3$W_Q15 <- df3$Q7*df3$Q15A_7

df3$Anx_Q15 <- df3$Anxiety_M*df3$Q15A_7


#Individual country OLS (Table 5)

#HongKong Set up

df5 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, HQ35, Q36, Q37, Q40, Q15A_7, xHQ41, Q44, W_Q15, Anx_Q15, Q31_HK, HCOUNTRY)%>%
  mutate(across(c("Q36", "Q44", "xHQ41", "Q31_HK"), factor))%>%
  filter(HCOUNTRY == 1)

#HongKongWorry (Model 5)

hk_wor <- lm(Q17_3 ~ Q7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 +  relevel(xHQ41, ref = "2") + Q31_HK, data= df5 )
summary(hk_wor)

kable(broom::tidy(hk_wor), digits = 2, align="c",
      caption="A data sample")

length(hk_wor$residuals)


#HongKongWorry interaction model (Model 15)

hk_wor_int <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +relevel(xHQ41, ref = "2") + W_Q15 + Q31_HK, data= df5 )
summary(hk_wor_int)

kable(broom::tidy(hk_wor_int), digits = 2, align="c",
      caption="A data sample")


#HongKongAnxiety (Model 6)

hk_anx <- lm(Q17_3 ~ Anxiety_M + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 + relevel(xHQ41, ref = "2") + Q31_HK, data= df5)
summary(hk_anx)

kable(broom::tidy(hk_anx), digits = 2, align="c",
      caption="A data sample")

length(hk_anx$residuals)


#HongKongAnxiety interaction model (Model 16)

hk_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + Anx_Q15 + Q31_HK, data= df5 )
summary(hk_anx_int)

kable(broom::tidy(hk_anx_int), digits = 2, align="c",
      caption="A data sample")


#Japan set up   

df6 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, HQ35, Q36, Q37, Q40, Q15A_7, xHQ41, Q44, W_Q15, Anx_Q15, Q31_JP, HCOUNTRY)%>%
  mutate(across(c("Q36", "Q44", "xHQ41", "Q31_JP"), factor))%>%
  filter(HCOUNTRY == 2)

#Japan Worry (Model 7)

jp_wor <- lm(Q17_3 ~ Q7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 + relevel(xHQ41, ref = "2") + Q31_JP, data= df6 )
summary(jp_wor)

length(jp_wor$residuals)

kable(broom::tidy(jp_wor), digits = 2, align="c",
      caption="A data sample")



#interaction model Japan Worry (Model 17)

jp_wor_int <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + W_Q15 + Q31_JP, data= df6 )

summary(jp_wor_int)

kable(broom::tidy(jp_wor_int), digits = 2, align="c",
      caption="A data sample")

#JapanAnxiety (Model 8)

jp_anx <-lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + Q31_JP, data= df6 )
summary(jp_anx)

length(jp_anx$residuals)

kable(broom::tidy(jp_anx), digits = 2, align="c",
      caption="A data sample")



#interaction model Japan Anxiety  (Model 18)

jp_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + Anx_Q15 + Q31_JP, data= df6 )
summary(jp_anx_int)

kable(broom::tidy(jp_anx_int), digits = 2, align="c",
      caption="A data sample")



#Singapore set up 

df7 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, HQ35, Q36, Q37, Q40, Q15A_7, xHQ41, Q44, W_Q15, Anx_Q15, Q31_SG, HCOUNTRY)%>%
  mutate(across(c("Q36", "Q44", "xHQ41", "Q31_SG"), factor))%>%
  filter(HCOUNTRY == 3)

#SingaporeWorry (Model 9)

sgp_wor <- lm(Q17_3 ~ Q7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 + relevel(xHQ41, ref = "2") + Q31_SG , data= df7  )
summary(sgp_wor)

length(sgp_wor$residuals)

kable(broom::tidy(sgp_wor), digits = 2, align="c",
      caption="A data sample")



#Singapore Worry interaction model (Model 19)

sgp_wor_int <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + W_Q15 + Q31_SG, data= df7  )

summary(sgp_wor_int)

kable(broom::tidy(sgp_wor_int), digits = 2, align="c",
      caption="A data sample")


#SingporeAnxiety (Model 10)

sgp_anx <-lm(Q17_3 ~ Anxiety_M  + Q15A_7  + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + Q31_SG ,data= df7  )
summary(sgp_anx)

length(sgp_anx$residuals)

kable(broom::tidy(sgp_anx), digits = 2, align="c",
      caption="A data sample")



#Singpore Anx interaction  model (Model 20)
sgp_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + Anx_Q15 + Q31_SG, data= df7  )
summary(sgp_anx_int)

kable(broom::tidy(sgp_anx_int), digits = 2, align="c",
      caption="A data sample")


#SouthKorea set up  

df8 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, HQ35, Q36, Q37, Q40, Q15A_7, xHQ41, Q44 ,W_Q15, Anx_Q15, Q31_SK, HCOUNTRY)%>%
  mutate(across(c("Q36", "Q44", "xHQ41", 'Q31_SK'), factor))%>%
  filter(HCOUNTRY == 4)

#SouthKorea Worry (Model 11)

sk_wor <- lm(Q17_3 ~ Q7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 + relevel(xHQ41, ref = "2") + Q31_SK, data= df8  )
summary(sk_wor)

kable(broom::tidy(sk_wor), digits = 2, align="c",
      caption="A data sample")



#South Korea Worry interaction model (Model 21)

sk_wor_int <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + W_Q15 + Q31_SK, data= df8  )
summary(sk_wor_int)

kable(broom::tidy(sk_wor_int), digits = 2, align="c",
      caption="A data sample")

#SouthKoreaAnx (Model 12)

sk_anx <-lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + Q31_SK, data = df8  )
summary(sk_anx)

length(sk_anx$residuals)

kable(broom::tidy(sk_anx), digits = 2, align="c",
      caption="A data sample")



#South korea Anx interaction  model (Model 22)
sk_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + Anx_Q15 + Q31_SK, data= df8  )
summary(sk_anx_int)

kable(broom::tidy(sk_anx_int), digits = 2, align="c",
      caption="A data sample")


#Taiwan set up   

df9 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, HQ35, Q36, Q37, Q40, xHQ41, Q15A_7, Q44, W_Q15, Anx_Q15, Q31_TW, HCOUNTRY)%>%
  mutate(across(c("Q36", "Q44", "xHQ41", "Q31_TW"), factor))%>%
  filter(HCOUNTRY == 5)

#Taiwan Worry (Model 13)

tw_wor <- lm(Q17_3 ~ Q7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + Q15A_7 + relevel(xHQ41, ref = "2") + Q31_TW, data = df9 )
summary(tw_wor)

kable(broom::tidy(tw_wor), digits = 2, align="c",
      caption="A data sample")



#Taiwan interaction model Worry (Model 23)

tw_wor_int <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + W_Q15 + Q31_TW , data = df9 )
summary(tw_wor_int)

kable(broom::tidy(tw_wor_int), digits = 2, align="c",
      caption="A data sample")

#Taiwan Anxiety (Model 14)

tw_anx <-lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2")  + relevel(xHQ41, ref = "2") + Q31_TW, data= df9 )
summary(tw_anx)

length(tw_anx$residuals)

kable(broom::tidy(tw_anx), digits=2, align="c",
      caption="A data sample")



#interaction Taiwan Anx model (Model 24)

tw_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + Anx_Q15 + Q31_TW, data= df9 )
summary(tw_anx_int)

kable(broom::tidy(tw_anx_int), digits=2, align="c",
      caption="A data sample")


##Create dummy variables 

df3$hk <- ifelse(df3$HCOUNTRY==1,  1, 0)

df3$jap <- ifelse(df3$HCOUNTRY==2,  1, 0)

df3$sing <- ifelse(df3$HCOUNTRY==3,  1, 0)

df3$kor <- ifelse(df3$HCOUNTRY==4,  1, 0)

df3$taiwan <- ifelse(df3$HCOUNTRY==5,  1, 0)

df3$TrustGov <- df3$Q16_1


## Fixed OLS Analysis (Table 4) 

df199 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M, Q15A_7, Anx_Q15, TrustGov, HQ35, Q36, Q37, Q40, xHQ41, Q44, HCOUNTRY, hk, jap, sing, kor, taiwan, W_Q15, Q16_1)%>%
  mutate(across(c("Q36", "Q44", "xHQ41"), factor))


#pooled anx (Model 3)

pooled_anx <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + TrustGov + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + jap + sing + kor + taiwan,  data = df199)
summary(pooled_anx)
kable(broom::tidy(pooled_anx), digits=2, align="c",
      caption="A data sample")

#pooled anx with Interaction effects (Model 4)

pooled_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + Anx_Q15 + TrustGov + HQ35 + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") + relevel(xHQ41, ref = "2") + jap + sing + kor + taiwan,  data = df199)
summary(pooled_anx_int)
kable(broom::tidy(pooled_anx_int), digits=2, align="c",
      caption="A data sample")


##pooled wor (Model 1)
pooled_wor <- lm(Q17_3 ~ Q7 + Q15A_7 + HQ35 + TrustGov + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2")  + relevel(xHQ41, ref = "2") + jap + sing + kor + taiwan, data = df199)

str(pooled_wor)

##pooled wor with interaction effects (Model 2)
pooled_wor_int<- lm(Q17_3 ~ Q7 + Q15A_7 + Q7*Q15A_7 + HQ35 + TrustGov + relevel(Q36,ref = '1') + Q37 + Q40 + relevel(Q44, ref = "2") +  relevel(xHQ41, ref = "2") + jap + sing + kor + taiwan, data = df199)
summary(pooled_wor_int)
kable(broom::tidy(pooled_wor_int2), digits=2, align="c",
      caption="A data sample")

str(pooled_wor_int)

#sensitivity analysis (via sensemakr) models 1-4 

pooled_wor_sens <- sensemakr(model = pooled_wor, treatment = "Q7", benchmark_covariates = "TrustGov", kd = 1:3) 

pooled_wor_int_sens <- sensemakr(model= pooled_wor_int2, treatment = "W_Q15", benchmark_covariates = "Q16_1", kd= 1:3)

pooled_anx_sens <- sensemakr(model = pooled_anx, treatment = "Anxiety_M", benchmark_covariates = "TrustGov", kd = 1:3)

pooled_anx_int_sens <- sensemakr(model= pooled_anx_int, treatment = "Anx_Q15", benchmark_covariates = "TrustGov", kd= 1:3)

plot(pooled_wor_sens) 
plot(pooled_wor_int_sens)
plot(pooled_anx_sens)
plot(pooled_anx_int_sens)

ovb_minimal_reporting(pooled_wor_sens, format= "html")
ovb_minimal_reporting(pooled_wor_int_sens, format= "html")
ovb_minimal_reporting(pooled_anx_sens, format= "html")
ovb_minimal_reporting(pooled_anx_int_sens, format= "html")



# Package interplot :: plot the effects of variables in Interaction term (Q7*Q15A_7 and Anxiety_M*Q15A_7 )

pooled_wor_int<- lm(Q17_3 ~ Q7 + Q15A_7 + Q7*Q15A_7 + HQ35 + Q16_1 + Q36 + Q37 + Q40 + Q44 + xHQ41 + jap + sing + kor + taiwan, data = df199)

interplot(m = pooled_wor_int, var1 = "Q7" , var2 = "Q15A_7" ) + ggplot2::labs(y = "Estimated Coefficient for Worry", x = "Percieved Border Stringency")  + ggtitle("Estimated coefficient of Worry and Perceived Border Stringency") + theme(panel.background = element_blank(), panel.grid.major = element_line(colour = "grey90"))

interactions::interact_plot(pooled_wor_int, pred = Q7, modx = Q15A_7, data = df199, modx.values = "plus-minus", legend.main = "Perceived border stringency") + labs(x = "Worry") + labs (y = "Support for immigration") + theme(legend.position = "bottom", axis.title = element_text(color = "black"), plot.title = element_text(color = "black"))


pooled_anx_int <- lm(Q17_3 ~ Anxiety_M + Q15A_7 + Anxiety_M*Q15A_7 + TrustGov + HQ35 + Q36 + Q37 + Q40 + Q44 + xHQ41 + jap + sing + kor + taiwan,  data = df199)

interplot(m = pooled_anx_int, var1 = "Anxiety_M" , var2 = "Q15A_7" ) + ggplot2::labs(y = "Estimated Coefficient for Anxiety", x = "Percieved Border Stringency") + theme(panel.background = element_blank(), panel.grid.major = element_line(colour = "grey90"))

interactions::interact_plot(pooled_anx_int, pred = Anxiety_M, modx = Q15A_7, data = df199, modx.values = "plus-minus" , legend.main = "Perceived border stringency") + labs( y = "Support for immigration") + labs(x = "Anxiety")+ theme(legend.position = "bottom", axis.title = element_text(color = "black"), plot.title = element_text(color = "black"))







##Casual Mediation Analysis, package :: mediation (Table 6)
#------1st mediation ( DV:17_3, IV: Q7, Mediator, Q6)


trace(mediation:::print.summary.mediate, 
      at = 11,
      tracer = quote({
        printCoefmat <- function(x, digits) {
          p <- x[, 4] 
          x[, 1:3] <- sprintf("%.6f", x[, 1:3])
          x[, 4] <- sprintf("%.2f", p)
          print(x, quote = FALSE, right = TRUE)
        } 
      }),
      print = FALSE)


df15 <- df3%>% 
  dplyr::select(Q17_3, Q7, Q6, Q5_4, Anxiety_M, Q16_1, HQ35, Q36, Q37, Q40, xHQ41, Q15A_7, Q44, hk, jap, sing, kor, taiwan)%>%
  mutate(across(c("Q36", "Q44", "xHQ41"), factor))

df16 <- na.omit(df15)

med.fit <- lm(Q6 ~ Q7 + HQ35 + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan, data= df16) 

out.fit <- glm(Q17_3~ Q7 + Q6 + HQ35 + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan,data= df16)

med.out <- mediation::mediate(med.fit, out.fit, boot = TRUE, treat = "Q7", mediator = "Q6", sims = 1000) 

summary(med.out)

mediation:::print.summary.mediate(summary(med.out))


#----sensitivity analysis for A 

sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 1000)

plot(sens.out, sens.par= "rho", main = "IV: Worry, Mediator: Economic threat")


#------B


med.fit2 <- lm(Q6 ~ Anxiety_M + HQ35 + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan, data= df16) 

out.fit2 <- glm(Q17_3~ Anxiety_M + Q6 + HQ35  + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan,data= df16) 

med.out2 <- mediation::mediate(med.fit2, out.fit2, boot = TRUE, treat = "Anxiety_M", mediator = "Q6", sims = 1000) 

summary(med.out2, sens.par= "rho", main = "Q6")

mediation:::print.summary.mediate(summary(med.out2))


#----sensitivity analysis for B

sens.out2 <- medsens(med.out2, rho.by = 0.1, effect.type = "indirect", sims = 1000)

plot(sens.out2, sens.par= "rho", main = "IV: Anxiety, Mediator: Economic threat")


#------C


med.fit3 <- lm(Q5_4 ~ Anxiety_M + HQ35  + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan, data= df16) 

out.fit3 <- glm(Q17_3~ Anxiety_M + Q5_4 + HQ35 + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan, data= df16)

med.out3 <- mediation::mediate(med.fit3, out.fit3, boot = TRUE, treat = "Anxiety_M", mediator = "Q5_4", sims = 1000) 

summary(med.out3)

mediation:::print.summary.mediate(summary(med.out3))


#----sensitivity analysis for C

sens.out3 <- medsens(med.out3, rho.by = 0.1, effect.type = "indirect", sims = 1000)

plot(sens.out3, main = "IV: Anxiety, Mediator: Impression towards foreigners")



#------D


med.fit4 <- lm(Q5_4 ~ Q7+ HQ35 + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan, data= df16) 

out.fit4 <- glm(Q17_3~ Q7 + Q5_4 + HQ35  + Q16_1 + Q36 + Q37 + Q40 + xHQ41 + Q15A_7 + Q44 + jap + sing + kor + taiwan,data= df16)

med.out4 <- mediation::mediate(med.fit4, out.fit4, boot = TRUE, treat = "Q7", mediator = "Q5_4", sims = 1000) 

summary(med.out4) 

mediation:::print.summary.mediate(summary(med.out4))


#----sensitivity analysis for D


sens.out4 <- medsens(med.out4, rho.by = 0.1, effect.type = "indirect", sims = 1000)

plot(sens.out4, sens.par= "rho", main = "IV: Worry, Mediator: Impression towards foreigners")


#Pooled OLS Ordinal (appendix)


df202 <- df3%>% 
  dplyr::select(Q17_3, Q7, Anxiety_M , HQ35, Q36, Q37, Q40, xHQ41, Q15A_7, Q44, HCOUNTRY, hk, jap, sing, kor, taiwan, W_Q15, Anx_Q15, Q16_1)

df202$Q17_3 = factor(df202$Q17_3, levels = c("1", "2", "3", "4", "5", "6", "7"), ordered = TRUE)
df202$Q36 = factor(df202$Q36, levels = c("1", "2"), ordered = TRUE)
df202$Q44 = factor(df202$Q44, levels = c("2", "1"), ordered = TRUE)
df202$xHQ41 = factor(df202$xHQ41, levels = c("2", "1", "3"), ordered = TRUE)
df202$HCOUNTRY = factor(df202$HCOUNTRY, levels = c("1", "2", "3", "4", "5"), ordered = TRUE)


pooled_wor_or_v2 <- polr(Q17_3 ~ Q7 + Q15A_7 + Q16_1 + HQ35 + Q36 + Q37 + Q40 + Q44 + xHQ41 + HCOUNTRY, data = df202, Hess=TRUE)
summary(pooled_wor_or_v2)

stargazer::stargazer(pooled_wor_or_v2, type="html", out="pooled_wor_or_v2.htm")


pooled_anx_or_v2 <- polr(Q17_3 ~ Anxiety_M + Q15A_7 + Q16_1 + HQ35 + Q36 + Q37 + Q40 + Q44 + xHQ41 + HCOUNTRY, data = df202, Hess=TRUE)
summary(pooled_anx_or_v2)
stargazer::stargazer(pooled_anx_or_v2, type="html", out="pooled_anx_or_v2.htm")



#Pooled OLS Ordinal analysis with Interaction term    

pooled_wor_int_or_v2 <- polr(Q17_3 ~ Q7 + Q15A_7 + Q16_1 + HQ35 + W_Q15 + Q36 + Q37 + Q40 + Q44 + xHQ41 + HCOUNTRY, data = df202, Hess=TRUE)
summary(pooled_wor_int_or_v2)
stargazer::stargazer(pooled_wor_int_or_v2, type="html", out="pooled_wor_int_or_v2.htm")


pooled_anx_int_or_2 <- polr(Q17_3 ~ Anxiety_M + Q15A_7 + Q16_1 + HQ35 + Anx_Q15 + Q36 + Q37 + Q40 + Q44 + xHQ41 + HCOUNTRY, data = df202, Hess=TRUE)
summary(pooled_anx_int_or_2)
stargazer::stargazer(pooled_anx_int_or_2, type="html", out="pooled_anx_int_or_2.htm")


##Correlation Matrix (Table 3) 

cor_matrix <- df3 %>%
  dplyr::select(Q17_3, Q7, Anxiety_M, Q15A_7, Q16_1, HQ35, Q36, Q37, Q40, Q44, xHQ41, Q5_4, Q6) 

head(cor_matrix, 6)

comp_cor <- cor(cor_matrix, method = "pearson", use = "complete.obs")

round(comp_cor, 4)


comp_cor_2 <- rcorr(as.matrix(cor_matrix))

comp_cor_2

flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
  )
}


flattenCorrMatrix(comp_cor_2$r,comp_cor_2$P)



